% =========================================================================
% run_simulation.m - (修改版) NMPC能量管理仿真主程序
%
% 核心修改:
% 1. 基于新的状态和控制向量，更新了初始条件和历史记录数组。
% 2. 在绘图部分增加了燃料电池温度和三个新的控制输入（流量）的可视化。
% 3. 【DEBUG 修正】修正了 lastMV 的初始值，以提供一个良好的初始猜测给求解器。
% =========================================================================

%% 1. 初始化
InitParameters; % 运行初始化脚本，创建所有参数和nlobj控制器

%% 2. 设置初始条件和扰动序列
% 初始状态: [SOC=80%, T_stack=353.15K, alpha=0, P_fc=min]
x0 = [0.8; 353.15; 0; FCParams.P_fc_min]; 
% 初始控制输入: [q_H2_in_min, q_air_min, q_W_min] (对应最小功率)
u0 = [FCParams.q_H2_min; FCParams.q_air_min; FCParams.q_W_min];

% 创建一个模拟的负载功率序列 (P_load) - 保持不变
N_sim = floor(SimParams.T_sim / SimParams.Ts); % 总仿真步数
P_load_sequence = 150 * ones(N_sim, 1);
P_load_sequence(50:150) = 400; % 5-15s: 高功率需求
P_load_sequence(151:175) = linspace(400,200,25);
P_load_sequence(176:300) = 200;% 15-30s: 中等功率
P_load_sequence(301:310) = linspace(200,450,10);
P_load_sequence(310:400) = 450;% 31-40s: 再次高功率
P_load_sequence(401:410) = linspace(450,150,10);
P_load_sequence(411:470) = 150; % 40-50s: 低功率
P_load_sequence(471:500) = linspace(150,0,30);
noise_sequence = 1+0.1*randn(size(P_load_sequence)); % 减小噪声影响
P_load_sequence = P_load_sequence.*noise_sequence;

%% 3. 准备仿真历史记录
x_history = zeros(N_sim+1, nx);
u_history = zeros(N_sim, nu);
x_history(1,:) = x0';
xk = x0; % 当前状态
lastMV = u0; % 【DEBUG 修正】为第一步提供合理的初始猜测

%% 4. 仿真主循环
disp('开始仿真... (程序将全速运行，请稍候)');

nloptions = nlmpcmoveopt; % 创建NMPC求解器选项
tic;
for k = 1:N_sim
    % --- a. 获取当前扰动 ---
    P_load_k = P_load_sequence(k);
    
    % --- b. 【核心】调用NMPC求解器计算最优控制量 ---
    % 传递P_load作为在线参数
    nloptions.Parameters = {P_load_k};
    % 获取上一步的控制输入，用于计算功率变化率
    [uk, nloptions, info] = nlmpcmove(nlobj, xk, lastMV, [], [], nloptions);
    
    % --- c. 将最优控制量应用到被控对象 ---
    % 注意：plant_update现在也使用详细模型
    xk_next = prediction_model(xk, uk, P_load_k);
    
    % --- d. 保存历史数据并更新状态 ---
    u_history(k,:) = uk';
    x_history(k+1,:) = xk_next';
    xk = xk_next;
    lastMV = uk; % 更新上一步控制量
end
sim_duration = toc;

disp('仿真结束。');
disp(['总耗时: ', num2str(sim_duration), ' 秒。']);

%% 5. 结果可视化 (已更新)
t = 0:SimParams.Ts:SimParams.T_sim;

% 绘制状态量
figure('Name', '状态量变化曲线');
subplot(4,1,1);
plot(t, x_history(:,1), 'b', 'LineWidth', 2);
title('SOC 变化曲线'); grid on; ylabel('SOC');
ylim([0.15, 0.85]);

subplot(4,1,2);
plot(t, x_history(:,2), 'Color', [0.8500 0.3250 0.0980], 'LineWidth', 2);
title('燃料电池温度'); grid on; ylabel('T_{stack} (K)');

subplot(4,1,3);
plot(t, x_history(:,4), 'r', 'LineWidth', 2); % P_fc现在是第4个状态
title('燃料电池输出功率'); grid on; ylabel('P_{fc} (W)');

subplot(4,1,4);
plot(t, x_history(:,3), 'g', 'LineWidth', 2);
title('燃料电池累积衰减'); grid on; ylabel('\alpha_{decay} (%)');
xlabel('时间 (s)');

% 绘制控制量
figure('Name', '控制量变化曲线');
subplot(3,1,1);
plot(t(1:end-1), u_history(:,1), 'm', 'LineWidth', 2);
title('控制量: 氢气流量'); grid on; ylabel('q_{H2} (mol/s)');

subplot(3,1,2);
plot(t(1:end-1), u_history(:,2), 'c', 'LineWidth', 2);
title('控制量: 空气流量'); grid on; ylabel('q_{air} (mol/s)');

subplot(3,1,3);
plot(t(1:end-1), u_history(:,3), 'k', 'LineWidth', 2);
title('控制量: 冷却水流量'); grid on; ylabel('q_{W} (mol/s)');
xlabel('时间 (s)');

% 绘制功率分配
figure('Name', '功率分配');
plot(t(1:end-1), P_load_sequence, 'k--', 'LineWidth', 2);
hold on;
P_fc_actual = x_history(2:end, 4); % 从状态历史记录中获取实际FC功率
P_comp = FCParams.P_comp_coeff * u_history(:,2).^2; % 简化的压气机功耗估算
P_bat_out = P_load_sequence - (P_fc_actual - P_comp);
area(t(1:end-1), [P_fc_actual, P_bat_out]);
title('功率分配'); grid on; ylabel('功率 (W)');
legend('总负载 P_{load}', '燃料电池 P_{fc}', '锂电池 P_{bat}');
xlabel('时间 (s)');
